Perturbative analysis of the conductivity in disordered mono— and bilayer graphene 
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FIG. 1: Two-loop corrections of the two-particle Green's 
function. Each diagram appears twice in the expansion. Solid 
lines in the upper/lower bow correspond to Green's functions 
of clean graphene {u — 0) G{ — z)/G{+z) and slashed lines to 
the disorder. 
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The DC conductivity of monolayer and bilayer graphene is studied perturbatively for different 
types of disorder. In the case of monolayer, an exact cancellation of logarithmic divergences occurs 
for all disorder types. The total conductivity correction for a random vector potential is zero, 
while for a random scalar potential and a random gap it acquires finite corrections. We identify 
the diagrams which are responsible for these corrections and extrapolate the finite contributions to 
higher orders which gives us general expressions for the conductivity of weakly disordered monolayer 
graphene. In the case of bilayer graphene, a cancellation of all contributions for all types of disorder 
takes place. Thus, the minimal conductivity of bilayer graphene turns out to be very robust against 
disorder. 
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Introduction. Monolayer graphene (MLG) represents 
a monoatomic sheet of carbon atoms arranged in a hon- 
eycomb lattice with lattice spacing a and next -neighbor 
hopping energy t w 2.8 eV. The transport properties of 
charge-neutral MLG are characterized by the semimetal- 
lic behavior with a point-like Fermi surface at two nodes 
(valleys) and linear low-energy dispersion in the vicinity 
of these valleys. This remarkable fact is one reason for 
the outstanding electronic properties of ML graphene 
0]. Perhaps the most prominent transport property of 
ML graphene is the minimal conductivity CTq = /hi: 
exactly at the Dirac point which has been observed in a 
number of experiments [H, [liB- Bilayer graphene (BLG) 
represents two ML honeycomb lattices with Bernal stack- 
ing, where the interlayer hopping processes are allowed 
with the energy tj^ w 0.4 eV. The main difference be- 
tween MLG and BLG is that the low-energy excitations 
of the latter have a quadratic spectrum in the vicinity of 
the valleys [5]. This difference causes a factor of two for 
the DC conductivity = 2e^//i7r. Experimentally, both 
values seem to depend only very weakly on disorder or 
thermal fluctuations [B-Qj which is supported by field- 
theoretical studies for various types of disorder l3" IB] ■ 

To the best of our knowledge, a systematic diagram- 
matic analysis of the conductivity of disordered graphene 
has not been performed so far. Usually, only certain 
types of diagrams are taken into account. Such approxi- 
mations cannot be considered as fully controllable, since 
each diagram in the perturbative expansion exhibits log- 
arithmic divergences. Therefore, the final result of cal- 
culations must crucially depend on a correct counting of 
diagrams, ft is the purpose of this paper to demonstrate 
this for the conductivity calculated within the Kubo for- 
malism. 

Kubo formula. Within the linear response theory, the 
conductivity of graphene per a spin and valley projection 
can be approximated for low frequencies (cj ^ 0) by the 
Kubo formula [l3| 




FIG. 2: Three-loop corrections of the two-particle Green's 
function. Each diagram appears four times. Generally, the 
combinatorial factor of each nth order diagram is 2". 
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is the analog of the (mean square) displacement func- 
tion of a classical random walk. The spin-valley de- 
generacy should be taken into account by multiplying 
it with an additional factor 4. The symbol r^. denotes 
the k component of the position operator; the brackets 
{■■■)g mean average with respect to disorder of strength 
g. The expression G^^,{z)Gj,,^{—z) is referred to as the 
two-particle Green's function and is commonly depicted 
diagrammatically as a closed loop <CZ> • Here, G denotes 
the one-particle Green's function with 



G^^{z) = hz + a - V + v'^Ur 



(3) 



where the vector a consists of Pauli matrices (Ti^2 and Uj. 
represents a disorder potential. The operator of kinetic 
energy reads Vj = ihvdi for MLG with the Fermi veloc- 
ity V = V^ta/2h, and Vi = fi^idf - a|)/2/i and V2 = 
h?^l^2/^J, for BLG, where ^ = 2t±h'^ /it^a^ denotes the 
electron band mass. For BLG we also define a character- 
istic velocity v = h/afj, = 3t'^a/2t±h which appears in the 
Green's function. Here, we consider three disorder types: 
1) random scalar potential, Ur = K-Cq! 2) random vector 
potential — a ■ and 3) random gap Ur = SMrCr^. 
We assume disorder to be Gaussian correlated with zero 
mean, i. e. (ii^) = 0, {ii^u^,) = {H/ v)'^ {g /mv'^)d{r — r') 
{{urUr') ~ {h/v)'^{g/mv'^)Sij6{r — r') for random vec- 
tor potential), with corresponding v and m = m^, the 
bare electron mass for MLG {m^v'^ « 6eV), and 
TO = /i for BLG (/zw^ ~ 30eV). Our assumption sug- 
gests g/{mv'^) ^ 1. Below we use a unit system with 



1, e^/h=l, 



TO„ 



1, 2/z = 1 and v = 1. 



Before we embark on the perturbative calculation, we 
briefly discuss the status quo of the field-theoretical ap- 
proach. In the case of a random vector potential, repre- 
senting random ripples in graphene, a bosonized replica 
approach yields for MLG a — ag The random gap 
case is of particular interest because it describes a metal- 
insulator transition due to the opening of the gap m. The 
displacement function Cg{m,ui) can be evaluated in this 
case by replacing the random field 5Mr by a more gen- 
eral random field. This mapping enables us to search for 
saddle points in a multidimensional manifold. It turns 
out that not just a single saddle point exists but a whole 
saddle-point manifold. As a result, we have one mass- 
less mode 1J|. The latter creates an singularity in 
Cg{m,ui), which cancels the factor in the conductiv- 
ity of Eq. ([J) . There are also massive modes around the 
saddle-point manifold. Altogether this leads to the fol- 
lowing scaling form; 

Cgim,io) = t±^!l)lc„{m/2,i, + 2irj)Kg, (4) 

CO 

with the displacement function of the pure system 



^ , ^ 2 1 

TT 4m^ - 
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Kg is the contribution of the massive modes. Although 
its form is unknown for g > 0, it is always finite with 
Kq = 1 for g — 0. The scattering rate r] is given by 
■q = {m? - ml)Q{ml - to^)/4 with rric = ke~'^l^ It 
should be noticed that the o;"^ singularity of Cg disap- 
pears for both MLG and BLG for 77 = and m > 0. In 
terms of the conductivity of Eq. ([T]) this yields for w ^ 
eventually 



—2 ) 



TO^) 
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This result indicates that a vanishes for any to > in 
perturbation theory because rric is always zero in the lat- 
ter. On the other hand, for to = the perturbation 
theory should reproduce a = a^Kg^ since the rric drops 
out of the conductivity in this case. 

Perturbative theory: Now we evaluate the displacement 
function perturbatively in powers of the disorder strength 
(7/2. In order to perform these calculations, the transla- 
tional invariance broken by disorder has to be restored by 
averaging over the Gaussian ensemble. Here, we perform 
this averaging using a diagrammatic representation. The 
first order conductivity corrections arise from the graphs 
depicted in Fig. [TJ Results of the evaluation of these di- 
agrams are shown in Table HI where we only retain con- 
stant contributions and logarithmically divergent terms. 
The second-order (three-loop order) diagrams contribut- 
ing to the conductivity are shown in Fig. [2] and results of 
their evaluation are summarized in Table [Til Technical 
details of the evaluation can be found in the supplemen- 
tary material 18|. The total combinatorial factor for 
each topological class of diagrams is 2". Each topologi- 
cal class of diagrams exhibits in turn further degeneracy 
due to diagram symmetries: The degeneracy factor of the 
first-order topological class 1 (Fig. [T]) and of the second- 
order topological classes 1, 2 and 7 (Fig. [2]) is one, while 
that of the first-order topological class 2 (Fig. [TJ and of 
the second-order topological classes 4, 5, 6 and 8 (Fig. [2]) 
is two and degeneracy factor of the second-order topolog- 
ical class 3 is four. All these factors should be carefully 
taken into account. 

First we discuss results for monolayer graphene. To the 
two-loop order, contributions to the displacement func- 



tion from all diagrams reveal divergences 



■InA/w. 



Here, A denotes the UV-cutoff and lo the imaginary fre- 
quency in contrast to Eqs. ^ and ([2|). The emergence 
of the logarithms is due to the divergence of the loop 
integrals 



<Pk A 
^k[z) ^ In- 



(2^)2 



(7) 



while the singularity uj ^ appears in the displacement 
function due to re-scaling of the factor in Eq. ([2]) [isj , 
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Diagram class Scalar disorder Gap disorder Vector disorder 

(m 1) X 1 a[l + 2^]/2 (2/3) -q[1 - 2^]/2 (2/3) 2c^^ (4/3) 

d 2) X 2 (-2/3) -at (-2/3) -2af (-4/3) 

Total X 2 a (0) -a (0) (0) 

TABLE I: Conductivity corrections from the first-order classes of diagrams (Fig. [T]) for MLG (BLG) in corresponding ctq 
units [l^. MLG results are shown for A ^ lj. Here we use the shorthands a — g/2-K, I — log(A/cj), and /3 = g/Sui. The 
degeneracy factors for each diagram and combinatorial factor 2 are taken into account (first column). 

Diagram class Scalar disorder Gap disorder Vector disorder 



(H 1) X 1 


a^[l + 4£ + 6f] 


/8 


(3/3^2) 


Q^[l -4^ + 6^^]/8 


(3/3V2) 






(6/3^) 


da 2) X 1 


IV2 




(5/3V6) 


5a^£Vl2 




(5/3V6) 






(10/JV3) 


(0 3) X 4 


-a'^£[l + e] 




(-3/3^) 


—a £ 




(-3/3^) 


-2a^£[l + 2^] 




(-12/3') 


(H 4) X 2 


a'^£[l~£]/2 




(2/3V3) 


q2^[1 -£]/2 




(2/3V3) 


2o?(\\ - I] 




(8/3V3) 


(0 5) X 2 






(0) 






(0) 






(0) 


(13 6) X 2 























(H 7) X 1 























(H 8) X 2 























Total X 4 


o?l2 




(0) 


o?l2 




(0) 










TABLE IL Conductivity corrections from the second-order classes of diagrams (Fig. [2]) for MLG (BLG) in corresponding oq 
units. All shorthands as above. The degeneracy factors and combinatorial factor 4 are taken into account (first column). 



but disappears in the conductivity because of the factor 
cj^ in Eq. ([T]), as in the field-theoretical approach dis- 
cussed above. At the three-loop level, contributions from 
diagrams containing intercrossing impurity lines (dia- 
grams 6, 7 and 8 in Fig. [5]) vanish after angular integra- 
tion. Contributions to the displacement function arising 
from each diagram with non-crossing impurity lines di- 
"^(In A/w)^, with some of them revealing sub- 



'InA/w. However, to both 



verge ^ 

dominant divergence 

orders g and g^, the sum over all conductivity contribu- 
tions is finite, i. e. singularities to both orders ~ InA/w 
and ^ (InA/oj)^ cancel each other exactly. For the ran- 
dom vector potential, the conductivity correction is zero 
to both two- and three-loop order. This is in accord with 
the findings of Refs. 
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Il5l . ll6l |. For random gap and ran- 
dom scalar potential to both orders, finite conductivity 
corrections are generated only by ladder diagrams, i. e. 
by diagrams 1 in Figs. [T] and [S] Provided that the can- 
cellation of singularities holds to higher orders as well. 



the analysis of higher order ladder diagrams yields the 
following general expression for the nth order {n ^ 1) 
conductivity correction: 



- Ml _ 
° 2"-i \2tt 



(8) 



with -t- for the random scalar potential and — for the 
random gap. This expression has been verified to the 
fourth order in perturbative expansion (five-loop order). 
The sum over n ^ 1 converges and we obtain for the 
conductivity 



>V,M 



1±A 

47r 

q- 



(9) 



for random scalar potential (+) and random gap (— ), cor- 
respondingly. It should be noticed that the conductivity 
is enhanced (reduced) by scalar-potential (gap) disorder. 



4 



This is plausible because the fluctuations of the former 
add particles to the system at the Dirac point, whereas 
the latter opens a fluctuating gap that should reduce the 
contribution to the conductivity. 

In the case of bilayer, the integral in Eq. ([7]) converges 
for A — > oo, giving / ^ const. Therefore, corrections 
to the displacement function arising from each nth order 
diagram become proportional to ui~'^~". Formally, this 



leads to the singularity 



of each nth order conduc- 



tivity correction for small w, as it is shown in Tables |T] 
and nil (values in brackets). Details of the evaluation of 
the first-order diagrams are summarized in [l8| . How- 
ever, the sum over all contributions gives a zero for all 
disorder types in both first and second order of perturba- 
tion theory. Provided this cancellation holds to higher or- 
ders as well, the minimal conductivity of bilayer graphene 
turns out to be very robust with respect to all types of 
disorder. 

Discussion: The conductivity formula Eq. ([9]) repre- 
sents the main result of our work. It has the required 
form a = &oKg. An important feature of this solu- 
tion is its scale invariance. Indeed, while this expression 
has been obtained for a finite frequency w after perform- 
ing the limit A oo, formally, the same result follows 
by keeping the cutoff A finite and performing the limit 
w — )■ 0. Hence, Eq. ([9]) can be regarded as an asymp- 
totically exact solution of the DC transport problem in 
disordered ML graphene. 

For the particular case of the random gap disorder in 
MLG with zero average gap, our results confirm find- 
ings of the recent numerical works [l^ H^l- For g ^ Q 
it reproduces the field theoretical result obtained in [l^] . 
The robustness of the minimal conductivity for BLG for 
zero average mass is in a remarkable agreement with the 
non-perturbative result obtained in [14| . On the other 
hand, the perturbation theory fails to describe a metal- 
insulator transition typical for a two-dimensional elec- 
tron gas [3, l20i] . This is because the scattering rate -q 
vanishes in perturbation theory, whereas a nonzero r] is 
the parameter which controls the metal-insulator transi- 
tion according to Eq. ^ . Hence, the area of applicability 
of perturbation theory is restricted to the metallic phase. 

Apart from the contribution to the displacement func- 
tion shown in Eq. ([2]) , there are two further contributions 
which emerge from the Kubo formula and might become 
important sufficiently far away from the Dirac point [l3| . 
In terms of one-particle Green's functions they corre- 
spond to the product G{±z)G{±z). However, close to the 
Dirac point these contributions can be neglected in com- 
parison to Eq. ([2]). Although the reason for neglecting 
them is not evident from the point of view of perturba- 
tion theory, the argument is provided by the field theory. 
The latter demonstrates the absence of massless modes 
at the saddle-point for these contributions [13: [13; i- ^■ 
at the Dirac point they are strongly suppressed and do 
not contribute to the transport. 



In conclusion, we have carefully studied perturbative 
corrections to the conductivity of disordered monolayer 
and bilayer graphene for different disorder types. Up to 
three-loop order we managed to show that in the case of 
ML graphene logarithmic divergences cancel each other 
exactly, irrespectively of the disorder type. Thus, the 
conductivity of weakly disordered monolayer graphene 
is modified by a finite correction. On the other hand, 
the minimal conductivity of bilayer graphene does not 
acquire any corrections for any type of disorder, as we 
have demonstrated for the first and second order in the 
perturbative expansion. 
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EVALUATION OF THE CONDUCTIVITY OF CLEAN MONOLAYER GRAPHENE 

Below we evaluate Kubo conductivity of clean monolayer graphene from 



where 



Cgico) ^r2Tr[G.o(^)Go.(-z)] = - ^ J] e-^^'^Tr{G.o(2)Go.(-z)} 



dpi 

with z — ie + uj/2. After performing Fourier transform we obtain 

CsW = -/(^ ^Tr{G,(z)Gg+p(-z)} 



p=0 



(1) 



(2) 



(3) 



p=0 



Now we perform Wick rotation z iz, (and lo ^ ico in the first term of Eq. ([I])), and obtain, after tracing over the 
spinor space, 



(Pq 9^ + g • (g + p) 



(27r)2 dpi [z2 + (q+p)2][z2 + g2] 



p=0 



The derivatives have to be performed first. We notice that 



and obtain 



and correspondingly 



^2 A{p) _ A"{p) A'{p)B'{p) A{p)[B'{p)f A{p)B"{p) 

' B3{p) B^{p) ' 

2^2 



dplB{p) B{p) B^ip) 
d"^ 2;2 + g • (g + p) 

9^ [z2 + (g+p)2][22 + ^2] 



p=0 



(z2 + g2)3 ■ 



27rz2 7rw2 



(4) 



(5) 



(6) 



(7) 



Inserting last expression into the Kubo formula Eq. ([T]) yields the minimal conductivity of ML graphene ctq = Xj'Ke^ jh. 
Evaluation for BLG goes analogously and yields CTq = Ij-ne^ jh. 

FIRST-ORDER PERTURBATIVE CONDUCTIVITY CORRECTIONS 

Below we calculate conductivity corrections from diagrams depicted in Fig. [TJ We start with the case of random 
gap for MLG. After Wick rotation, first diagram reads 



and the second: 



<Pq f ffk 92 

T^Tr {G^+fc(-iz)Gg+fc+p(iz)cr3Gg+p(iz)G^(-iz)cr3} 



2 J (27r)2 J (27r)2 dp, 
g [ <Pq f (fk 92 



p=0 



(2^)2 J (2^)2 dp\ 



■Tr{G^+p(iz)Gq(-iz)cr3Gq+fc(-iz)cr3Gg(-iz)} 



(8) 



(9) 



p=0 



2 




with Green's functions 



(10) 



Here we use the slashed notation ^ = a ■ q. We start with Eq. Using the fact that each integration runs from 
—oo to +00 we may shift k ^ k — q. Then the fc-integration acts only on the third Green's function and diverges 
logarithmically 



A2 



(11) 



where the UV-cutofF A is introduced only temporarily and will be removed at the end of calculation. This quantity 
couples to the unity matrix and can be pulled out from the trace operator. After performing the trace, the g-integral 
becomes 



d^q 



2iz 



1 



P 



(z2 + (5+p)2) 



" 4^ z3(A2 + z2)2 ' 



(27r)2 dpliz^ + q'f 
such that the contribution to the displacement function from this diagram becomes 

g A2(A2 + 2z2) 



167r2 z2(A2+22)S 



In 1 



A2 

^2 



(12) 



(13) 



The evaluation of the expression in Eq. ([S]) is more cumbersome. First, we can also reshift the momentum k ^ k — q 
as above. Performing some reordering in the numerator we arrive at the following expression: 



'-'9,1 



(1) _ 9 



(Pq f cPk Tr[z^ + z^{ff + ^^ + f) + ^i^+f)}+^^ + f){^+^)^] 



2 J (2^)2 y (2^)2 dpi {z' + k'){z' + q')(z' + {p + qy){z' + {p + ky) 



p=0 



Using contraction formula for Pauli-matrices 



the trace over the last term in the numerator becomes 

Trm^ + lf)(f + m = 2{q^k^ + q^k-p + k^q-p + p^k ■ q). 

Thus we get from Eq. (HI]) 

d2fc f (Pq d"^ z^ + z^ [p^ + k ■ {k+p) + q- {q+p)\+ q^k^ + q^k -p + k^q-p + p^k ■ q 



(2^)2 y (2^)2 dpi 



{z' + k'){z' + q'){z' + (p + qY){z' + (p+ fc)^) 



(14) 

(15) 
(16) 

(17) 



p=0 



Using Eq. ([5]) we obtain for each single term 
A'{p) _ 2(z2 + fc . g) 



B{p) 



A'{p)B'{p) 



B^ip) 



p=o {z' + q'nz' + k'f' B^p) 
_4Mz' + q^) + q,{z' + k')]\ A{p) 



{z' + k^iz' + qy ' B^ip) 



B"{p) 



^ 2Mz^ + q^)+q,{z^ + k')f 
{z' + kYiz' + qY ' 

_ 2[(z2+g2)+4fc^5^ + (^2^^2)] 



p=0 



(z2 + A;2)2(z2+g2)2 



3 



Now we are in the position to perform integrals: 



9 



(fq A"{p) 



(2^ 



-9 



-2 



(27r)2 J (27r)2 B{p) 
A'{p)B'{p) , A{p)[B'{p)f 



BHp) 
(fk r (Pq 



(27r)2 J (27r) 



S3(p) 

A{p)B"{p) 
BHp) 



p=0 



p=0 



87r2 z2(A2 + 22)2' 

g 

"87r2 z2(A2+z2)2 

A2 



In 1 



A2 



5 



47r2 z2(A2 + ^2) 



In 1 



A= 



The contribution to the displacement function from the first-order ladder diagram then reads: 



A^ 



9 A2(A2 + 2z2) / A2 

B7r2z2(A2 + z2)2 87r2 z2(A2 + z2)2 \ 



(18) 



Taking combinatorial factors into account, the total first-order perturbative contribution to the displacement function 
becomes 



•'3,2 



A4 



47r2 z2(A2 + z2)2' 



(19) 



where all terms containing logarithms eliminate each other. After removing the cutoff A — >■ cx3 we obtain a finite 
conductivity correction 



2 ^' 27r2 27r ' 



(20) 



For the case of random scalar disorder we have to replace Pauli matrices (Tg by unity matrices in Eqs. (|8l) and ([9]). 



This does not imply any changes on Cg ^ but changes {. At the level of Eqs. (|14p and p7)) this can be seen by the 
different sign of the term z^p^ or z^p^, correspondingly, while all other terms remain the same. Then it leads to the 
change in 



A"{p) 



B{p) 



2{k-q-z^) 



p=0 



(g2+z2)2(^2^^2)2' 



with all other terms remaining without any changes. Eventually, we will obtain a finite conductivity correction as in 
Eq. (|20|) but with the opposite sign. 

For the case of random vector disorder we have to replace Pauli matrices (Xg by cr^, fi — 1,2 in Eqs. (|8]) and ([9]) 
and perform an additional summation over /i. Let's evaluate Eq. ^ for this case. Shifting k ^ k — q we separate the 
singular part from the regular one. The singular part remains the same as in Eq. (jlll) . Evaluation of the regular part 
goes in line with the analogous calculation for the random gap disorder with the only difference that contraction of 
Pauli matrices gives an additional factor 2. Hence, the result for the conductivity correction reads 



,(2) _ 



27r2 I A* / w ' 



(21) 



where we redefined the cutoff by A ^ \/A^ 
has to notice that 



up- 12. For the evaluation of Eq. ([8|) for the case of vector disorder one 
Y^a^fa^ = 0, (22) 



2{t + 0. 



(23) 



This leads after performing the trace to the following expression: 



r^i^ - -a 

<^o,i - 9 



(Pq ^2 (z2 + g2)(2;2 + + (^2 + p)q . p {^z'^ q^)k ■ p - p^q ■ k + 2k ■ p q ■ p 



Pk 



(27r)2y (27r)2 dpf 



(^2 + q2)[z^ + P)(22 ^_ (p ^_ g)2)(22 + + ;,)2) 



(24) 



33=0 
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After taking derivatives with respect to the momentum pi we get 



A"{p) 



B(p) 
A{p)[B'{p)Y 



B^ip) 



_ Ak^q^ -2q-k 

~ (^2 + ^2)2(^2 + ^2)2; 



A'{p)B'{p) 



BHp) 
A{p)B"{p) 



= 2 



p=0 



p=0 



+ ^2)3(^2 _|_ p)3 



B^{p) 



p=0 



[g,(z2 + fc2)+fc,(^2+^2)]2 
(z2 + g2)3(22 + p)3 ^ 

(z2 + fc2) + 4g.fc. + (^2 + g,2) 
(z2 + g2)2(^2 + p)2 ' 



With exception of the first expression all of them are the same as compared to the case of random gap disorder. 
Hence, the character of the singularity does not change, while the regular part averages to zero after performing 
angular integration. With the same redefinition of the cutoff as above, the result for this conductivity contribution 
reads 



,(1) _ 9 



The total contribution of all diagrams to the conductivity then becomes 



2a 



(1) 



4a 



(2) 



0. 



"1 — ^'^l I 

Now we turn our attention to BLG. Green's function for BLG reads 



=F iz 



z^ + q'i 



(25) 



(26) 



(27) 



with Qi = ?! — q2 and Q2 = 2q^q2. For the gap disorder, the evaluation of Eq. © goes analogously to the MLG 
calculations. After reshifting of momenta, the ^-integration can be performed separately which gives 

d^k 



{2n) 



G,i-iz) 



(28) 



This factor can be pulled out of the trace operator and product of sigma matrices , corresponding to the disorder 
type, gives a unity matrix. Therefore, the result of the evaluation will be the same for the scalar disorder too, while 
that for the vector disorder will be twice this value. After performing the trace over the remaining terms we arrive at 



d^q 92 



i2z 



z^-q^ + 2[(gi + pi)2 - q^]{ql - qj) + 8{q^ + Pi)qiql 



(27r)2 9p2 (^2 + ^4)2 z-^ + [{qi+Pi 

Performing the derivatives with the help of Eq. ([5]) we obtain exploiting the rotational invariance 



Pi=0 



92 z^-q^ + 2[{q^ +J5i)2 - ql]{ql - gf) + 8(^1 + Pi)qiql 



dpi z"^ + [{qi+PiY + ql? 

The contribution to the displacement function from this diagram then reads 

d^q q2 g 



8g2 



Pi=0 



z2 + ^4 ■ 



C. 



(1) 

3,2 



(27r)2 (z2+g4)3 



27rw3' 



(29) 



(30) 



(31) 



for both gap and scalar disorder, and twice this value for the vector disorder. The conductivity correction for gap/scalar 
disorder is therefore proportional to lo^^: 



,(1) 



So; ' 



(32) 



where the minimal conductivity of BLG is CTq = 2/t:. For the vector disorder it must be multiplied by 2. 

The type of disorder matters if we evaluate the contribution to the displacement function from the first diagram 
in Fig. [T] given analytically in Eq. ([5]). In order to simplify the calculation we reshift momenta as in the MLG case 
discussed above. In terms of Eq. (O 



(1) _ 



dl^q f d^k 92 ^(p) 



(27r)2 J (27r)2 dpj B{p) 



(33) 



p=0 
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we obtain for all disorder types the same denominator 

B{p) = {z' + q'){z^ + k'){z^ + {q+pf){z^ + {k+pf), (34) 

but different expressions for numerators. For scalar/gap disorder it can be written as 

A{p) = 2{z'' + z''[Q{q)-Q{q+p)+Q{k)-Q{k+p) 

T {Q{q) ■ Qik) + Qiq+p) ■ Q{k+p) ~ Qiq) ■ Q{k + p) - Q{k) ■ Q{q + p))] 

+ Q{q) ■ Q{q + P) Q{k) ■Q{k+p)-Q{q)- Q{k + p) Q{k) ■ Q{q + p) + Q{q) ■ Q{k) Q{q + p) ■ Q{k + p)} 

with (-)/(+) for scalar/gap disorder. Using the relations cr^'^cr^ = and a ^Q{q)Q{q+p)(j ^ = 2Q{q+p)Q{q), 
the numerator for the vector disorder becomes after performing the trace 

A{p) = A{z^ + z''[Q{q)-Q{q+p) + Q{k)-Q{k+p)] 

+ Qiq) -Qiq+p) Qik) ■ Qik + p) - Qiq) ■ Qik) Qiq+p)- Qik + p) + Qiq) ■ Qik + p) Qik) -Qiq+p)}. 

After performing the derivatives, taking the limit p — > and carrying out integrations we obtain for scalar/gap 
disorder in terms of Eq. ([5]): 



(fiq [ (Pk A"ip) 



<Pq 

(27r)2 J (27r) 



(27r)2 J (27r)2 
(Pk fA'ip)B'ip) Aip)[B'ip)] 



0, 



B^ip) 



BHp) 



(Pq f (Pk Aip)B"(p) 



(27r)2 J (27r)2 B^{p) 



p=0 



p=0 



p=0 



1 



'i'KZ^ ' 

1 



27rz3 



Thefore, the contribution to the displacement function for scalar/gap disorder becomes with z — u /2 



and the corresponding conductivity correction 



In the case of vector disorder we have 



a, ^^jji 



,(1) _ 5C^0 



4iu 



(Pq f Pk A"(p) 



-2 



(27r)2 J (27r)2 Bip) 
Pq f Pk (A'ip)B'ip) Aip)[B'ip)f 



= 0, 



(27r)2 J (27r)2 



B2(p) B^ip) 
Pq f Pk Aip)B"ip) 



(27r)2 7 (27r)2 
which yields for the displacement function correction 



B^ip) 



p=0 



p=0 



1 



27rz3' 
1 



0.1 ~ ^ ■ 



and eventually for the conductivity correction 



2uj 



In all cases, the total first-order conductivity correction is zero: 



(1) 



4(7, 



(1) 



(35) 
(36) 
(37) 

(38) 

(39) 

(40) 
(41) 
(42) 

(43) 
(44) 
(45) 
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FIG. 2: Momentum reshifting of the second-order ladder diagram. 



EVALUATION ALGORITHM FOR HIGHER-ORDER DIAGRAMS 

The higher-order Feynman graphs were evaluated with the help of computer algebra package Mathematica. We 
used the algorithm described below: 

• Draw diagrams and distribute momentum indices in accord with momentum conservation; 

• Reshift momenta appropriately, perform the trace and introduce the polar coordinates; 

• Perform derivatives with respect to pi and send p — > 0. For higher-order diagrams it is convenient to use Eq. ([5]); 

• Perform angular integrals. At this stage all second-order diagrams containing crossed impurity lines vanish; 

• Perform radial integrals from to the upper cutoff A. Replace z = uj/2 and redefine for MLG the cutoff 
via A — >■ \/A^ — uP- 12. For BLG, the limit A — > oo is carried out. Use the command FullSimplify[ expression, 
Assumptions— > {A > 0, uj > 0}]; 

• Sum up contributions from all diagrams taking correct degeneracy factors into account, Use the command 
FullSimplify[e2:pression, Assumptions— >■ {A > 0, w > 0}]. At this stage all divergent contributions drop out; 

• Multiply with — a;^/2(g/2)", where n denotes the diagram order and perform the limit A — > oo. Only the finite 
corrections survive. 

Example 1: As a first example we consider the second-order ladder diagram for MLG with random scalar disorder 
shown in Fig. [2j 

1. First we define propagators. The Mathematica code reads: 

• Define Pauli matrices: 

= {{1,0}, {0,1}} 

a, = {{0,1}, {1,0}} 

= {{0,-1}, {1,0}} 

• Define propagators for positive and negative frequencies: 

GreenPl = (-1 zctq + ql (Ti + q2 o-2)/(z^ + ql^ + q22) 
GreenMi = ( I zcto + ql (Ti + q2 cr2)/(z^ + ql^ + q22) 

2. Now we draw the diagram and put momentum indices as required by momentum conservation. This is displayed 
by the left graph in Fig. [2j Next we reshift momenta in two steps: 1) t t — k and 2) k k~ q. The result 
can be seen on the right graph in Fig. [21 The Mathematica code which implements these steps reads: 

Bl = GreenMi 

B2 = GreenPl/. {ql ^ ql + pi, q2 ^ q2 + p2} 

B3 = GreenPl/. {ql ^ ql + kl + pi, q2 ^ q2 + k2 + p2} 

B4 = GreenPl/. {ql ^ ql + kl + tl + pi, q2 ^- q2 + k2 + t2 + p2} 

B5 = GreenMi/. {ql ^ ql + kl + tl, q2 ^ q2 + k2 + 12} 

B6 = GreenMi/. {ql ^ql + kl, q2 ^ q2 + k2} 
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reshift momenta p 




q+k+p 

FIG. 3: Momentum reshifting of the second-order X-diagram. 

Perform the trace, reshift momenta as discussed above, put p2 — > and introduce polar coordinates 

Al = SimpUfy [(((Dot[Bl, B2, B3, B4, B5, B6] [[1, 1]] + Dot[Bl, B2, B3, B4, B5, B6] [[2, 2]])/(2 Pi)^ 
/. {tl ^ tl - kl - ql, t2 t2 - k2 - q2}) 
/. {kl ^ kl - ql, k2 ^ k2 - q2}) 

/. {p2 ^0, ql ^ q Cos[x], q2 ^ q Sin[x], kl ^ k Cos[y], k2 ^ k Sin[y], tl t Cos[u], t2 t Sin[u]}] 

3. Perform derivatives with respect to pi, put pi — > and perform angular integrations 

A2 = Integrate[Integrate[Integrate[Simplify[D[Al,pl,pl]/.pl 0],{x, 0, 2Pi}],{y, 0, 2Pi}],{u, 0, 2 Pi}] 

4. Perform radial integrations from to A: 

A3 = (Integrate[t A2, t]/. t ^ A) - (Integrate[t A2, t]/. t 0) 
A4 = (Integrate [k A3, k]/. k ^ A) - (Integrate[k A3, k]/.k^O) 
A5 = (Integrate[q A4, q]/. q ^ A) - (Integratefq A4, q]/. q ^ 0) 



5. Simplify, multiply with — aj^/2, reshift the cutoff by A — > \/A^ — and simplify again 
A6 = FullSimphfy[A5, Assumptions ^ A > 0, z > 0] 

A7 = FullSimplify[-a;V2 (A6/. {z ^ uj/2, A ^ Sqrt[A2 ~ uj^]/2}), Assumptions ^ {A > 0, w > 0}] 
The result of the last operation reads 

(A2 - c^2)3 + 2A2 In - f 2(A2 - uj^f + SiA" - o.^) In - 

(jJ \ UJ 



87r3A6 

which upon multiplication with (5/2)^ and retaining only finite and divergent terms for A — > 00 yields the conductivity 
correction 

2 . or . n2N 



(1) / 1 r 1 , A 3 

="M8(2^+2(2^^"lJ + 4 



27r Lu 



where (Tq = I/tt, as it can be seen in Table II in the main text. 

Example 2: As another example we demonstrate vanishing of the so-called diagram with maximally crossed 
impurity lines shown in Fig. [3] for MLG. 

1. First we distribute momenta in accord with the momentum conservation as shown on the left in Fig. [3] and 
reshift them in the following steps: 1) t ^ t — k — q, 2) q ^ q — p. The result is shown on the right in Fig. [S] 
To implement this with Mathematica we have to define the propagators: 

Bl = GreenMi 

B2 = GreenPl/. {ql ^ ql +pl, q2 ^ q2 + p2} 

B3 ^ GreenPl/. {ql ^ ql +kl + pl, q2 ^ q2 + k2 + p2} 

B4 = GreenPl/. {ql ^ ql + kl + tl + pi, q2 ^ q2 + k2 + t2 + p2} 

B5 = GreenMi/. {ql ^ ql + kl + tl, q2 q2 + k2 + t2} 

B6 = GreenMi/. {ql ^ ql +tl, q2 ^ q2 + 12} 

multiply them, perform the trace and reshift momenta as shown above: 

B7 = Simplify [(((Simplify [Dot [Bl, B2, B3, B4, B5, B6] [[1, 1]] + Dot[Bl, B2, B3, B4, B5, B6] [[2, 2]]]) 
/. {tl -J> tl - kl - ql, t2 ^ t2 - k2 - q2})/. {ql ^ ql - pi, q2 -J> q2 - p2})] 
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2. Now we put p2 to zero and introduce polar coordinates: 

Tl = Simplifyp?/. {p2 ^0, ql ^ q Cos[x], q2 ^ q Sin[x], kl ^ k Cos[y], 
k2 ^ k Sin[y], tl t Cos[u], t2 t Sin[u]}] 

3. In order to perform differentiations with, respect to pi we use Eq. ([5]). Tlie numerator (A) and denominator (B) 
of Tl are separated out using following commands: 

A = Numerator [Tl] 
B = Denominator[Tl] 

4. Now we generate first term from Eq. ([5|), put pi — > and perform integration over the angle y from to 2tt: 

T2 ^ Simplify[(Integrate[Simplify[D[A, pi, pl]/B/. {pi -> 0}], y]/. y ^ 2Pi) 
- (Integrate[Simplify[D[A, pi, pl]/B/. {pi ^ 0}], y]/.y^O)] 

The result of this integration reads 

47rz^ cos[u + x] 



qt (q2 + z2)2(t2+z2)2 
and the next integration with respect to x or u yields zero. 

5. Generate second term from Eq. ([5]), put pi — > and perform integration over the angle y: 

T3 = Simplify[(Integrate[Simplify[D[A, pl]D[B, pl]/B2/. pi ^ 0], y]/.y->2Pi) 

- (Integrate[Simplify[D[A, pl]D[B, pl]/B2/. pi ^ 0], y]/.y^O)] 

The result is zero. 

6. Generate third term from Eq. ([5]), put pi — > and perform integration over the angle y: 

T4 = FullSimplify[(Integrate[Simplify[(A(D[B, pl])2/B3)/.pl 0], y]/.y^2Pi) 

- (Integrate[Simplify[(A(D[B, pl])2/B3)/.pl ^ 0], y]/. y 0)] 

The result is zero. 

7. Generate fourth term from Eq. ([5]), put pi — >■ and perform integration over the angle y: 

T5 = Simplify[(Integrate[Simplifv[(A/B2) D[B, pi, pl]/.pl-^0], y]/.y-^2Pi) 

- (Integrate[Simplify[(A/B2) D[B, pi, pl]/.pl^O], y]/.y^O)] 

The result is zero. 



